\documentclass[UTF8]{ctexart}

\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amsthm}
\usepackage{xeCJK}
\usepackage{graphicx} %插入图片的宏包
\usepackage{float} %设置图片浮动位置的宏包
\usepackage{subfigure} %插入多图时用子图显示的宏包
\usepackage{geometry}
\usepackage{enumitem}
\usepackage{verbatimbox}
\usepackage[ruled,linesnumbered]{algorithm2e}
\usepackage[colorlinks,linkcolor=blue,citecolor=blue]{hyperref}
\usepackage{cite}

\newtheorem{example}{例}             % 整体编号
\newtheorem{theorem}{定理}[section]  % 按 section 编号
\newtheorem{definition}{定义}
\newtheorem{axiom}{公理}
\newtheorem{property}{性质}
\newtheorem{proposition}{命题}
\newtheorem{lemma}{引理}
\newtheorem{corollary}{推论}
\newtheorem{remark}{注解}
\newtheorem{condition}{条件}
\newtheorem{conclusion}{结论}
\newtheorem{assumption}{假设}

\geometry{a4paper,scale=0.8}
\CTEXsetup[format={\Large\bfseries}]{section}

\title{\textbf{Manderbrot Set 的生成和探索}}

\date{}

\author{\CJKfamily{kai} 黄文\hbox{\scalebox{0.6}[1]{羽}\kern-.1em\scalebox{0.5}[1]{中}}}

\begin{document}

\maketitle

\begin{abstract}
    \CJKfamily{kai} 本文借助自编的带deflate无损压缩的PNG编码器, 实现了Mandelbrot集指定区域图像的生成, 以某一点为中心, 逐渐缩小半径, 生成了一系列Mandelbrot集局部图像. 最后介绍了本文采用的PNG编码器的原理, 并分析了压缩效果.
\end{abstract}

\section{引言}

本文的目的是生成一系列Mandelbrot集图像. 第一部分介绍了Mandelbrot集的生成原理, 并给出了一种染色方法, 给出了算法流程图描述, 并在Github仓库中提供了算法实现的c++源代码. 

同时, 本文作者注意到无压缩图像格式对存储空间的巨大浪费, 因此研究了PNG图像的编码格式, 借助Github用户 \href{https://Github.com/queensun}{queensun} 用C语言实现的zlib压缩算法, 用c++完成了一个实现简单、封装优雅的PNG编码器. 本文的后半部分着重介绍了PNG的一种最简单的支持无损压缩的编码格式, 并介绍了本文实现的带deflate无损压缩的PNG编码器, 同时也在Github仓库中提供了编码器的c++源代码. 

\section{问题背景}

Mandelbrot集是复动力系统中著名的图形, 由Mandelbrot教授在20世纪70年代发现, 被称为“上帝的指纹”. 其特点是无限分形. 本文作者第一次接触到Mandelbrot集是在尹永成教授的分析学课上. 尹老师展示了他在复动力系统中建设性的研究成果, 并向学生展示了Mandelbrot集的美丽图形. 

Mandelbrot集的绘制一直是计算数学家和计算机学家所热衷的事情, 由于Mandelbrot集本身的无限分形性质, 只要染色函数足够优美, 就可以通过局部放大生成一系列美观的图形, 这些图形可以广泛应用于设计领域. 

\section{Mandelbrot集的原理与生成}

\subsection{Mandelbrot集的数学原理}

\begin{definition}
    记$f_c(z)=z^2+c$, 其中$c\in\mathbb{C}$. 记$f_c^{\circ(0)}(z)=f_c(z)$, 迭代地定义$f^{\circ(k)}_c(z)=f_c(f^{\circ(k-1)}_c(z))\;(k\geq 1)$. 称集合
    \begin{equation*}
        \mathcal{M}=\{c\in\mathbb{C}:f_c^{\circ(n)}(0)\text{关于}n\text{有界}\}
    \end{equation*}
    为\textbf{Mandelbrot集}.
\end{definition}

\begin{property}
    设$c\in\mathbb{C}$, 若$c\in\mathcal{M}$, 则$|f_c^{\circ(n)}(0)|\leq 2\;(\forall n\in\mathbb{N})$.
\end{property}

证明见参考文献\cite{branner1989mandelbrot}. 这个定理可以推导出Mandelbrot集的生成算法. 

\begin{definition}
    对于$c\notin \mathcal{M}$, 令
    \begin{equation}
        L(c)=\min\{n:f_c^{\circ(n)}(0)|> 2\}
    \end{equation}
    称$L(c)$为点$c$的\textbf{层数}, 记第$n(n\geq 0)$层的所有点构成的集合为$L_n$, 即: 
    \begin{equation}
        L_n=\{c:L(c)=n\}
    \end{equation}
\end{definition}

\begin{property}
    $\forall n\geq 1$, $L_n$是一个连通集, 其连通数为2, 外边界为$\{c:|P_{n}(c)|=2\}$, 其中$P_n$是一个$n$次多项式.
\end{property}

性质2表明, $L_n$具有简单光滑的封闭边界. 其具体形式见参考文献\cite{张国栋2005Mandelbrot}. 由性质2易得以下结果.

\begin{property}
    $\forall c\notin \mathcal{M}$, 存在一个充分小的邻域$N(c)$, 使得在该邻域内, 函数$L(z)$要么恒为常值, 要么只有两种取值且由一条简单光滑曲线分割两个取值不同的区域.
\end{property}

性质3表明, 任何Mandelbrot集外部的点, 都只能“有限分形”. 这两条性质将指导Mandelbrot集局部放大的绘图中心点选取.

\begin{property}
    设$\{z_n\}$是一列复平面上的点列, 且$z_n\in L_n$, 则其任何子列都存在聚点, 且聚点$z^*\in \partial \mathcal{M}$.
\end{property}

\subsection{Mandelbrot集的生成算法}

利用性质1, 可以得知, 对于一个点$c$, 令$z_0=0,\;z_n=z_{n-1}^2+c$, 不断地执行迭代, 当$|z_n|>2$时即可判定$c\notin \mathcal{M}$. 然而计算机不可能达到无穷次迭代, 因此设置迭代次数上限$N$, 当$|z_n|\leq 2\;(\forall n=0,...,N)$时, 就近似地认为$z_n\in \mathcal{M}$. 显然$N$越大计算结果越精确, 但计算代价也越大. 

对于$c\notin \mathcal{M}$, 给出着色函数$g:\mathbb{N} \to \mathbb{Z}_{255}^3$, 将$c$点的颜色设为$g(L(c))$, 这样即可获得分层的彩色图像. 本文采取的着色函数为: $g(n)=(255\log_2(n)\mod 256,\;155\log_2(n)\mod 256,\;150)$

下面给出局部染色算法, 算法中的\textbf{PixelIndexToCodinate}是像素点坐标到轴坐标的转换函数: 

\begin{algorithm}[H]
	\caption{drawMandelbrotImage}
	\small
	\KwIn{ 中心点坐标center, 绘图区域直径diam, 图片宽度W, 图片高度H, 最大迭代次数N, 着色函数g, 文件名Fname}
	\KwOut{ picture.png }
	
	\For{$i=1$ to $H$ }
	{
        \For{$j=1$ to $W$}
        {
            $c$ = \textbf{PixelIndexToCodinate}(center,diam,$i,j$)\\
            $z_0$ = 0\\
            $k$ = 0\\
            \While{$|z_k|\leq 2$ and $k<N$}
            {
                $z_{k+1}=z_k^2+c$\\
                $k\gets k+1$\\
            }
            \eIf{$k$ == $N$}
            {
                setPixel($i$,$j$,(0,0,0))
            }
            {
                setPixel($i$,$j$,$g(k)$)
            }
        }
	}
    \textbf{print image to} Fname.png
\end{algorithm}

\subsection{Mandelbrot集的局部放大算法}

Mandelbrot集最为人熟知的性质就是无限分形, 以特定的点为中心, 不断放大, 可以生成无穷无尽的图像. 下面给出一个算法, 对于给定的绘图中心$c$, 生成一共28张图像, 每张图都是前一张图放大2.73倍的结果. 

\begin{algorithm}[H]
	\caption{getLocalImages}
	\small
	\KwIn{ 中心点坐标center, 图片宽度W, 图片高度H, 着色函数g}
	\KwOut{ picture\_1.png, ..., picture\_32.png }
	
    $d=5.0$\\
    $N=250$\\
	\For{$i=1$ to $28$ }
	{
        \textbf{drawMandelbrotImage}(center,d,W,H,N,g,picture\_$i$)
        $d\gets d \times 0.99^100$
        $N\gets round(N \times 1.0025^100)$
	}
\end{algorithm}

上述算法最大的困难在于中心的选取, 中心选取不当产生的图会非常单调. 

根据性质3, 对于一个Mandelbrot集外部的点, 它总是有限分形的, 放大到一定倍数后整个图将只具有一种颜色. 而对于Mandelbrot集内部的点, 放大到一定倍数后整个图中所有点都位于Mandelbrot集内部, 也不再具有分形性质. 因此, 最理想的绘图中心是Mandelbrot集的边界点. 

Mandelbrot集的边界点无法用初等表达式给出, 在计算机中无法精确表示. 因此根据性质4, 一个好的策略是, 选择一个浮点数$c$, 使得$L(c)$充分大, 从而充分趋近于Mandelbrot集边界, 进而获得优美的图案. 

\subsection{Mandelbrot集的图像展示}

通过局部放大、手动调整中心、再局部放大、再手动调整中心……不断重复, 得到一个充分好的点：$c_0=(-0.72624686941319,0.240376999706965)$

以$c_0$为中心, 执行局部放大算法(Algorithm.2), 获得下面一系列Mandelbrot集局部图像. （附件中的视频展示了放大过程）

\begin{figure}[H]
    \centering
    \includegraphics[width=0.18\textwidth]{img/pic0.png}
    \includegraphics[width=0.18\textwidth]{img/pic50.png}
    \includegraphics[width=0.18\textwidth]{img/pic100.png}
    \includegraphics[width=0.18\textwidth]{img/pic150.png}
    \\
    \centering
    \vspace{3pt} \includegraphics[width=0.18\textwidth]{img/pic200.png}
    \includegraphics[width=0.18\textwidth]{img/pic250.png}
    \includegraphics[width=0.18\textwidth]{img/pic300.png}
    \includegraphics[width=0.18\textwidth]{img/pic350.png}
    \\
    \centering
    \vspace{3pt} \includegraphics[width=0.18\textwidth]{img/pic400.png}
    \includegraphics[width=0.18\textwidth]{img/pic450.png}
    \includegraphics[width=0.18\textwidth]{img/pic500.png}
    \includegraphics[width=0.18\textwidth]{img/pic550.png}
    \\
    \centering
    \vspace{3pt} \includegraphics[width=0.18\textwidth]{img/pic600.png}
    \includegraphics[width=0.18\textwidth]{img/pic650.png}
    \includegraphics[width=0.18\textwidth]{img/pic700.png}
    \includegraphics[width=0.18\textwidth]{img/pic750.png}
    \\
    \centering
    \vspace{3pt} \includegraphics[width=0.18\textwidth]{img/pic800.png}
    \includegraphics[width=0.18\textwidth]{img/pic850.png}
    \includegraphics[width=0.18\textwidth]{img/pic900.png}
    \includegraphics[width=0.18\textwidth]{img/pic950.png}
    \\
\end{figure}

\begin{figure}[H]
    \centering
    \vspace{3pt} \includegraphics[width=0.18\textwidth]{img/pic1000.png}
    \includegraphics[width=0.18\textwidth]{img/pic1050.png}
    \includegraphics[width=0.18\textwidth]{img/pic1100.png}
    \includegraphics[width=0.18\textwidth]{img/pic1150.png}
    \\
    \centering
    \vspace{3pt} \includegraphics[width=0.18\textwidth]{img/pic1200.png}
    \includegraphics[width=0.18\textwidth]{img/pic1250.png}
    \includegraphics[width=0.18\textwidth]{img/pic1300.png}
    \includegraphics[width=0.18\textwidth]{img/pic1350.png}
    \caption{以(-0.72624686941319,0.240376999706965)为中心展开的Mandelbrot集图像}
\end{figure}

\section{无损压缩PNG编码器的原理与实现}

\subsection{PNG编码方式}

本文实现的是一个简单的8位真彩色图像带压缩PNG编码器. 为简明起见, 本文采用符合需求的最简单编码方式: 由文件头和IHDR, IDAT, IEND三个数据块组成. 其中文件头为固定8字节: 

\begin{verbatim}
    8950 4e47 0d0a 1a0a
\end{verbatim}

一个数据块有长度（Length）、标签（Tag）、数据（Data）、校验码（CRC）四个部分, 说明如下.

\begin{table}[H]
    \centering
    \begin{tabular}{|c|c|c|}
        \hline
        \textbf{名称} & \textbf{长度} & \textbf{说明} \\ \hline
        Length & 4 bit & Data部分的字节数, 用32位无符号整型方式编码 \\ \hline
        Tag & 4 bit & 数据块名称对应的ASCII码的十六进制, 例如IHDR的Tag部分为 \verb |4948 4452| \\ \hline
        Data & 可变长 & 数据块实际存储的数据内容 \\ \hline
        CRC & 4 bit & 从Tag开头到Data结束部分的4字节CRC校验码 \\ \hline
    \end{tabular}
\end{table}

IEND部分不存储数据, 其Data部分为空, 因此其CRC校验码也恒为常值. 下面说明IHDR与IDAT的Data部分编码.

\subsubsection{IHDR的Data编码}

IHDR的Data部分由13个字节构成, 其具体含义如下: 

\begin{enumerate}[itemindent=2em]
    \setlength{\itemsep}{-5pt}
    \item \verb| width |: 4bit, 32位无符号整数表示图片宽度
    \item \verb| height |: 4bit, 32位无符号整数表示图片高度
    \item \verb| bit depth |: 1bit, 表示图像深度, 本文取 \verb |0x08| 表示8位
    \item \verb| color type |: 1bit, 表示色彩类型, 本为取 \verb |0x02| 表示真彩色图像
    \item \verb| compression |: 1bit, 表示压缩方式, 本文取 \verb |0x00| 表示deflate/inflate无损压缩
    \item \verb| filter |: 1bit, 表示滤波器方法, 本文取 \verb |0x00| 表示自适应选择, 实际上本文不做滤波器
    \item \verb| interlance |: 1bit, 表示扫描方法, 本文取 \verb |0x00| 表示非隔行扫描
\end{enumerate}

\subsubsection{IDAT的Data编码}

对于一个$w\times h$大小的图像, 其IDAT的Data部分原始数据有$h*(1+w)$ bit, 设$R_{ij}$, $G_{ij}$, $B_{ij}$分别表示第$i$行第$j$列像素点的R/G/B值, 其编码方式可表述为: 

\newpage
\begin{equation*}
    0,R_{11},G_{11},B_{11},R_{12},G_{12},B_{12},\cdots,R_{1w},G_{1w},B_{1w},0,R_{21},G_{21},B_{21},R_{22},G_{22},B_{22},\cdots,R_{2w},G_{2w},B_{2w},\cdots
\end{equation*}

将编码存入 \verb |originData| 中, 然后调用zlib开源库函数: 

\begin{verbatim}
    compress(compressData, &compressLen, originData, h*(1+w))
\end{verbatim}

获得压缩后的数据, 存储在 \verb |compressData| 中, 其长度为 \verb |compressLen|, 压缩后的数据构成IDAT的Data部分的全部内容.

\subsection{PNG编码器的实现与接口}

本文根据上述原理实现了带deflate无损压缩的PNG编码器, 源码见: \href{https://Github.com/EbolaEmperor/MathSoftware/tree/main/Homework-4}{Github - 源码链接}, 函数接口见 \verb |png.h| , 函数实现见 \verb |png.cpp|. 使用示例参考绘制Mandelbrot集的主程序. 

\subsection{测试结果}

图1中的所有图像均由本文PNG编码器生成.其中每个图都是$1920\times 1080$尺寸的高清图像, 若不带压缩, 每个图的大小约为: $\frac{1920\times 1080\times 3}{1024\times 1024}=6.2$(MB).使用本文PNG编码器, 28个图像总计52.0 MB, 平均单个图片大小1.86MB, 平均压缩率为$70\%$.

事实上, 由于压缩算法的普遍特性, 当图像的重复元素较多时, 压缩率较高, 反之较低.图1中压缩率最大与压缩率最小的两个图像如下. 左图仅152.3 kB, 压缩率为$97.6\%$；右图2.6MB, 压缩率为$58.1\%$.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.2\textwidth]{img/pic0.png}
    \hspace{15mm}
    \includegraphics[width=0.2\textwidth]{img/pic750.png}
    \caption{图1中的两个典型图像}
\end{figure}

\section{总结}

本文实现了Mandelbrot集的绘制, 并得出了一系列优美的图像. 这些图像不仅能帮助复动力系统的理论研究者进行研究, 也能在工业设计、纺织等领域得到广泛的应用. 此外, 本文的代码实现优雅, 易读易用, 也不失为一个亮点. 

美中不足之处是, 由于本文的代码实现中采用double类型的浮点数运算, 只能精确到小数点后大约18位, 因此无法继续将图形放大下去, 本文几乎已经在double的精度范围内放大到了极限. 这个问题可能需要引入高精度才能解决. 

\bibliographystyle{plain}
\bibliography{reference}

\end{document}